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Introduction Riboswitches are cis-acting regulatory RNA elements prevalently located in the 
leader sequences of bacterial mRNA. An adenine sensing riboswitch cis-regulates adeninosine deami¬ 
nase gene {add) in Vibrio vulnificus. The structural mechanism regulating its conformational changes 
upon ligand binding mostly remains to be elucidated. In this open framework it has been suggested 
that the ligand stabilizes the interaction of the distal “kissing loop” complex. Using accurate full- 
atom molecular dynamics with explicit solvent in combination with enhanced sampling techniques 
and advanced analysis methods it could be possible to provide a more detailed perspective on the 
formation of these tertiary contacts. 

Methods In this work, we used umbrella sampling simulations to study the thermodynamics of 
the kissing loop complex in the presence and in the absence of the cognate ligand. We enforced 
the breaking/formation of the loop-loop interaction restraining the distance between the two loops. 

We also assessed the convergence of the results by using two alternative initialization protocols. A 
structural analysis was performed using a novel approach to analyze base contacts. 

Results Contacts between the two loops were progressively lost when larger inter-loop distances 
were enforced. Inter-loop Watson-Crick contacts survived at larger separation when compared with 
non-canonical pairing and stacking interactions. Intra-loop stacking contacts remained formed upon 
loop undocking. Our simulations qualitatively indicated that the ligand could stabilize the kissing 
loop complex. We also compared with previously published simulation studies. 

Discussion and Conclusions Kissing complex stabilization given by the ligand was compatible 
with available experimental data. However, the dependence of its value on the initialization protocol 
of the umbrella sampling simulations posed some questions on the quantitative interpretation of the 
results and called for better converged enhanced sampling simulations. 


Introduction 

Riboswitches are portions of ribonucleic acid (RNA) able to regulate gene expression in bacteria and plants at several 
levels. They bind their sensed ligands without the need for protein factors. To regulate their target gene, riboswitch 
can either act on transcription, on translation, or, more rarely, as interfering, antisense or self-splicing RNAs [T]. More 
precisely, riboswitches are cis-acting RNA elements prevalently located in the leader sequences of bacterial mRNA [2] 
that regulate the expression of the same gene from which they have been transcribed. They are composed of an aptamer 
domain that binds the effector ligand and of an expression platform, usually located downstream of the aptamer, that 
transduces the ligand-induced conformational switch into the gene expression regulation [a 11]. Riboswitches are 
classified according to the nature of the sensed ligand [1]. Among them, the purine-sensing riboswitches emerge as 
important model systems for exploring various aspects of RNA structure and function [5] because of their structural 
simplicity and relatively small size. Within the purine family the add adenine-sensing riboswitch (A-riboswitch) is 
one of the most characterized. Found in the mRNA 5’-untraslated region, it cis-regulates the adenosine deaminase 
gene in Vibrio vulnificus acting rho-independently at the translational level [6]. Its regulatory activity depends on 
the availability of the ligand: in the presence of adenine the riboswitch is in the ON state, and the protein synthesis 
is permitted, whereas in the absence of the ligand the riboswitch folds into the OFF state blocking the translation 
initiation (Figure 1). The ligand-bound structure of its aptamer [71|8] is a junction of three stems (PI, P2, P3) with 
the ligand completely encapsulated into the structure (Figure 2). There are three structurally important regions: the 
binding pocket, the Pl-stem and the loop-loop tertiary interaction between L2 and L3, usually called “kissing loops”. 
The latter includes two inter-loop Watson-Crick (WC) base pairs [9l[T0]. The ligand-dependent structural mechanism 
inducing the switch between the ON- and the OFF-state in the A-riboswitch mostly remains to be elucidated. The role 
of the ligand in the structural organization of the aptamer has been investigated using structure-based fluorescence 
spectroscopy El, multidimensional NMR techniques m and single-molecule experiments m- These investigations 
however lack both the atomistic details and the distinct energetic contributions associated to ligand binding. In this 
open framework, in particular, it has been suggested experimentally that the ligand stabilizes the interaction of the 
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FIG. 1: Mechanism of action of A-riboswitch Cartoon showing the OFF (left) and ON (right) states of the A-riboswitch. 
When ligand is not present the ribosome binding site (orange, RBS) is paired with a portion of the aptamer and translation is 
blocked. When ligand is present the RBS is free to interact with the ribosome and translation can be initiated. 



FIG. 2: Structure of the add aptamer domain A) Three dimensional representation of the aptamer with the adenine 
bound. The stems are shown in grey and labeled. The backbone of the loops and of the junctions is shown in orange. B) 
Close-up on the loop-loop (L2 and L3) interaction with focus on the Watson-Crick base pairs (G25-C49, G26-C48, in blue). C) 
Close-up on the ligand binding site with the adenine (red) paired with U62. 


distal kissing complex M- At the same time a stable kissing interaction seems to contribute to the ligand binding 
energy stabilizing the complex in a cooperative fashion miiiiiE]. Also in silico techniques have been used obtaining 
an accurate description of the system from a structural point of view HMH]. 

In a few cases a computational approach has been employed to provide a thermodynamic characterization of the 
system HUH]. In particular, Allner et al. [20] computed the free-energy profile corresponding to the formation of 
the kissing complex using molecular dynamics (MD) simulations in explicit solvent, both in the presence and in the 
absence of the ligand, using the CHARMM36 force field [22[|23]. MD does not require experimental inputs and can 
in principle be used in a predictive fashion. However, accuracy of atomistic force fields is still debated and it is thus 
very important to compare results obtained employing different sets of parameters. 

In this paper we use atomistic MD with the latest variant of the Amber force field m in combination with enhanced 
sampling techniques [25] to provide a more detailed perspective on the formation of the kissing loop complex. The 
combined approach allows this contribution to be dissected from the other ligand-aptamer interactions and the impact 
of the ligand on the stability of the loop-loop interaction to be quantified. We reproduce exactly the same protocol that 
has been used by Allner et al. m in order to perform a fair comparison between the two force-fields on this particular 
system. Effects of the initialization protocol on the results of umbrella sampling simulations are also discussed in 
detail. 
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Methods 

In this work, we used umbrella sampling (US) simulations [26] to study the thermodynamics of the kissing loop in 
the presence and in the absence of the cognate ligand. We enforced the breaking/formation of the loop-loop interaction 
steering the distance between the two loops and then used the resulting structures as starting conformations for US 
with multiple restraints m- Simulations were carried out with the Gromacs 4.6.3 program package [28| combined 
with the PLUMED 2.0 plug-in [29]. All the simulation parameters are discussed in detail in the following subsections. 


System set up 

All simulations reported hereafter were performed on two systems: the add aptamer domain complexed with adenine 
(Holo form) and without adenine bound (Apo form). In both cases we used the X-ray structure solved by Serganov 
et al. [PDB:1Y26] [ 2 ]. The ligand was removed to simulate the unbound state. MD simulations were performed using 
the Amber99 force field [30] refined with the parmbscO a /7 corrections [31] and the latest y torsional parameters [24] . 
The general Amber force field [32] was used to parametrize the ligand. Partial atomic charges were assigned using the 
restricted electrostatic potential fit method [33] based on an electronic structure calculation at the HF/6-31G* level of 
theory performed with Gaussian03 [34] . The electrostatic interactions were calculated using the particle-mesh Ewald 
method [35] and bond-lengths were constrained with LINGS [36]. The systems were set-up following exactly the 
protocol described in Allner et al. [ 20 ] : aptamers were solvated in a rhombic dodecahedron having 8 nm as box vector 
lenghts, with a Mg^+-H 2 0 solution using approximately 11000 TIP3 molecules [37], and a recent parametrization 
for divalent cations [38]. The 5 crystallographic Mg^+ were initially kept at their respective position, whereas the 
additional 30 ions added to neutralize the system ([Mg^+] = 0.18 M) were randomly placed. A steepest descent 
minimization (150 steps) was performed followed by 200 ps of MD at constant temperature (298 K, using stochastic 
velocity rescaling [39]) and pressure (1 atm, using the Berendsen barostat [40]) with positional restraints on both 
RNA and ions so as to equilibrate water. This procedure was repeated first removing the constraints on the ions 
and then removing all the remaining constraints. Finally, 12 ns unrestrained simulations at constant volume were 
performed for each system. 


Umbrella sampling 

In order to compute the thermodynamic stability of the loop-loop interaction we employed US simulations with 
multiple harmonic restraints. The distance between the center of mass (CoM) of the backbone atoms of the two 
loops (L 2 : bases 20-26; L3: bases 48-54; Figure 2 B) was used as a collective variable (CV). We will refer hereafter 
to this distance as L. 44 uniformly spaced reference values were taken in the range spanning from 12.5 to 34 A, 
and restraints with stiffness /c = 20 {kcal/mol)/A? were employed. In the production phase of the US simulations 
each of the 44 windows was run for 5 ns. A very important issue in US simulation is the generation of the starting 
conformations. We here performed two independent US sampling simulations, using starting conformations generated 
with two different protocols (hereafter referred to as forward and backward). To generate the starting points for 
the forward US simulations we employed the same protocol used by Allner et al. m- Namely, we ran a series of 
44 short (0.25 ns) simulations with a stiffer restraint {k = 40 {kcal/mol)/K^) keeping the CV at the 44 reference 
values, where each simulation was initialized from the last frame of the previous one. In this way, before each US 
window starting structure was sampled, we let the system equilibrate. The reference values were iterated allowing an 
increasing distance between the loops. The backward US simulation was initialized with an equivalent procedure but 
iterating the restraints in the opposite order, i.e. starting from the structure with undocked loops. In principle, if US 
simulations are converged, the result should be independent of the initialization procedure. 


Analysis methods 

The data were analyzed using the last 4 ns of each window. The potential of mean force (PMF) profiles were 
constructed using the weighted histogram analysis method (WHAM) [27] implemented by Grossfield [41] taking the 
CV values distribution resulting from the US simulations. This implementation of WHAM allows to compute errors 
with a bootstrapping procedure that assumes uncorrelated samples. To avoid artifacts due to possible correlations we 
instead adopted a blocking procedure. Namely, we split the final 4 ns of each trajectory in four blocks of 1 ns each 
and performed the WHAM calculation using only a single block from each simulations. The four resulting profiles 
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FIG. 3: Potential of mean force for kissing-complex formation Potential of mean force (PMF) as a function of the 
distance between the centers of mass of the L2 and L3 loops. Results for Holo and Apo forms are shown as obtained from 
two independent umbrella sampling simulations using different protocols to obtain the initial structures (forward, Fwd, and 
backward, Bwd, see main text for dehnition). Fwd and Bwd profiles are aligned at the maximum distance. 


are aligned at their CV starting value (12.1 A for the forward profiles, 34.4 A for the backward profiles) and error at 
each point is computed as the standard deviation among the four profiles divided by \/4. 

To define the number of stacking interactions and the number of base-pair contacts a local coordinate system 
was constructed in the centre of each six-membered rings, with the x axis pointing towards the C2 atom and the z 
axis orthogonal to the ring plane. The pairing and stacking relationship between two bases j and k is based on the 
vector Vjk, i.e. the position of ring center k relative to the coordinate system constructed on base j. The criteria for 
determining the canonical WC base pairs are: 1) the base pair must be AU or GC; 2) The relative position of the 
bases is compatible with the geometry of a WC interaction. The latter condition is considered satisfied when the 
product of the Gaussian function M{rjk] x M{rkj] > 10“^. Mean /i and covariance cf were obtained from 
the empirical distribution of WC pairs in the crystal structure of the large ribosomal subunit [42]. The criteria for 

determining the non-canonical base pairs are: 1) the ellipsoidal distance Vj^ = ^ 

and Vkj < a/2 ^; 2) \zjk\ and \zkj\ < 2 A; 3) it is not a WC pair. The criteria for determining the stacking base 
pairs are: 1) the ellipsoidal distance Vj^ < and Vj^j < 2) \zjk\ and \zkj\ > 2 A; 3) x‘^-j^ + < 5 A and 

x\j + yf.- < 5 A. This procedure yields similar results compared to the MC-annotate software [43] and was shown to 
be useful for characterizing both structural and dynamical properties of RNA molecules |44|. The software used to 
perform this structural analysis is available online (http://github.com/srnas/barnaba). 


Results 

Forward process 

The analysis of the Holo forward and Apo forward US trajectories allowed the PMF for the disruption of the kissing 
complex to be computed. The resulting profiles are plotted in Figure 3 for both Holo and Apo systems. The PMF 
shows a minimum at L 12.5 A, corresponding to the initial structure. The free energy change upon disruption 
of the kissing complex for the Holo structure is AG = 52 ± 2 kcal/mol. For the Apo structure the stability of the 
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FIG. 4: Analysis of inter and intra-loop interactions Average count of inter and intra-loop interactions from umbrella 
sampling simulations. Results are shown for both Apo and Holo forms, using both forward (Fwd) and backward (Bwd) protocols 
(see main text for definition). Watson-Crick and non-Wat son-Crick pairings as well as intra and inter-loop stackings are shown 
as indicated. 


complex is largely reduced to AG = 35 ± 3 kcal/mol. The stabilization of the kissing complex provided by the ligand 
can thus be estimated as AAG = 17 ± 3 kcal/mol. To understand which are the interactions that are relevant for 
the kissing complex formation we analyze inter-loop pairings and inter- and intra-loop stacking interactions for each 
of the restrained simulations (Figure 4). 

For both the Holo and the Apo forms, at a L ^ 16 A, only the two inter-loop WC base pairs (G25-C49, G26-C48) 
peculiar of the loop-loop interaction are still formed. On the contrary, all the non-canonical base pairs are disrupted. 
In the Apo structure the inter-loop WC pairings were irreversibly lost at L > 23 A, whereas in the Holo structure they 
are at least partially maintained until L ^ 30 A. We also analyzed the rupture of stacking interactions, distinguishing 
intra-loop and inter-loop contacts. Inter-loop stacking behaves in a manner qualitatively similar to the inter-loop 
WC pairings, going to zero at a distance L « 23 A (Apo) and 30 A (Holo). On the contrary, the intra-loop stacking 
interactions are still present when the kissing loop is disrupted, indicating that the internal structure of the two loops 
is preserved during undocking. It can be observed that in the Holo simulation the number of intra-strand stacking 
slightly decreases (~ 5) for 29 A < L < 30 A because of the distortion in the structure induced by the one of the two 
inter-loop WC pairings. After this residual interaction is lost all the intra-loop stacking contacts are recovered. 


Backward process 

In order to better assess the convergence of the free energy landscape for kissing complex formation, we also 
reconstructed the PMF profiles of the Holo and Apo structure from US simulations initialized with the backward 
process (Figure 2). The forward and backward profiles were aligned at L = 34 A, since at that distance the starting 
structure of the backward process is equal to the final structure of the forward one. The free-energy change upon 
docking is estimated taking the difference between the minimum value of the PMF (Holo: L ^ 19 A; Apo: L ^ 16 A) 
and its value for the undocked structure (L = 34 A): for the Holo AG = -3.2 ± 0.9 kcal/mol, for the Apo AG = -5.9 
± 0.6 kcal/mol. Albeit negative, these numbers are too small and not compatible with the ones found in the forward 
process. This is a clear signature of hysteresis in the pulling procedure that strongly biases the initial starting points 
of the US simulation. The reason for this discrepancy can be better understood by performing a structural analysis of 
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the interactions on the different US windows. As it can be seen in Figure 4, in the backward process the native WC 
base pairs are not reformed. In general, a few contacts are formed between the two loops but they are not enough to 
stabilize the kissing complex. To be sure that this is a systematic effect we also tried a few alternative settings for 
backward simulations. Results are presented in the Appendix. 


Discussion and Conclusions 

Our calculations provide quantitative and atomistic details on the mechanism of kissing loop breaking and formation 
in the add riboswitch aptamer domain. The results can be compared with those recently obtained by Allner et al 
[20] on the same system using the CHARMM36 force field [22l [23] . In particular the free-energy computed with the 
forward process has been obtained with an identical protocol so as to allow a fair comparison between the force fields. 
In our work the estimated stability of the kissing loop complex is AG = 52 ± 2 kcal/mol (Holo) and AG = 35 ± 
3 kcal/mol (Apo), so that upon ligand binding AAG=17 kcal/mol > 0. On the contrary, Allner et al. reported 
AAG = —10 kcal/mol < 0. The sign of AAG indicates whether the ligand binding and the formation of the kissing 
loop complex are cooperative (positive) or anticooperative (negative). Results obtained with the two force fields thus 
interestingly lead to two opposite pictures. 

Recent experiments probed the differential ligand affinity in aptamers with mutations hindering the formation of 
the kissing complex M- The change in affinity indicates a cooperativity between ligand binding and kissing complex 
formation. This stabilization has been estimated to be AAG ~ 6 kcal/mol. This number should be interpreted with 
caution since it is based on the assumption that the mutated aptamer mimics a ligand-bound state that is accessible 
to the wild type aptamer M- Results obtained with Amber force field are in qualitative agreement with this picture. 

Recently, the thermodynamics of other stand-alone kissing complexes has been also studied using different bio¬ 
physical techniques [13 Sg. In these two experimental works the stability of the loop-loop interactions resulted to 
be in the range 8 < AG < 14 kcal/mol. Stability depends on the exact sequence and set of intra-loop interactions, 
but is always on the order of ten kcal/mol. The estimated stability of the kissing loop complex in our calculation, 
namely AG = 52 ± 2 kcal/mol and AG = 35 ± 3 kcal/mol for Holo and Apo respectively, is thus much larger than 
expected. Results obtained with CHARMM indicate a lower AG for both the systems m, in better agreement with 
experimental results, even if still overestimated. Our result could be affected by the known overestimation of stacking 
interactions in the Amber force field [47j. Additionally, we would like to point out that the overestimation found 
in our calculations could also be a consequence of difficult convergence in the US simulations. To test if the US 
simulation are effectively converged, we tried to recover the profiles from simulations initialized with the backward 
process, with a procedure inspired by two directional pulling in steered MD [48H5Q] . The discrepancy between forward 
and backward process is an index of high dependence of the PMF on the initialization procedure and poses some 
questions on the actual convergence of the US simulations. Similarly to steered MD, one can expect that simulations in 
the forward and backward process are respectively overestimating and underestimating the kissing complex stability. 
Optimal results can be obtained in steered MD by combining simulations performed with both protocols IHISO]. We 
stress here that even if the forward simulations apparently recover the qualitative behavior of the general accepted 
model, they cannot be trusted for a quantitative estimation of the free-energy change. The fact that the backward 
process cannot reach the native docked state is a signature of a barrier in an orthogonal degree of freedom that is not 
properly sampled. A possible candidate is the barrier related to the desolvation of the loops, required to form the 
correct interstrand interactions. Additionally, we observe that pulling on the distance between the two loops does not 
necessarily induce the entropic reduction required upon docking. These issues are expected to affect both forward 
and backward pulling. Our simulation could not give an estimate of the additional barriers, but we can assume that 
these issues equally affects the Holo and the Apo systems. Thus, the converged AAG upon ligand binding should 
be somewhere in between results from the forward and backward simulations. Thus we can expect the AAG to be 
in a wide range between -2.7 kcal/mol and 17 kcal/mol, which is in qualitative agreement with already mentioned 
experiments m- 

The convergence problem is not related to the US method itself but to the difficulty of describing such a complex 
docking event using a single distance as a CV. This variable is not sufficient to drive the system through the appropriate 
transition states. This is likely due to the existence of additional barriers on hidden degrees of freedom (e.g. solvation). 
We believe that in order to reliably quantify the AAG for this system with US or other biased sampling methods one 
should employ more complex CVs which are closer to the actual reaction coordinate. 

In conclusion, in this work we addressed the formation of the kissing loop complex in the A-riboswitch aptamer by 
means of accurate molecular dynamics simulations in explicit solvent combined with enhanced sampling techniques 
in presence or absence of the cognate ligand. Results are compatible with experiments and suggest that the ligand 
stabilizes the kissing loop formation. However, our results also spot some weakness of the umbrella sampling method 
and call for calculations performed with more advanced techniques, which will be the subject of future investigations. 
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FIG. 5: Analysis of interactions with alternative backward protocol Count of inter and intra-loop interactions in the 
44 starting snapshots resulted from the different Apo US simulation initialization procedures. Results are shown for the 4 Apo 
runs, using the backward and forward protocols, both with k = 40{kcal/mol)/ (I, II, respectively in orange and pale orange), 
and the two alternative backward protocols with the softer restraint (III in blue from protocol A, and IV in light blue from 
protocol B). Watson-Crick and non-Watson-Crick pairings as well as intra and inter-loop stackings are shown in the different 
panels. 


Appendix 

In order to assess the backward procedure for the US method, we repeated it for the Apo form using a softer 
restraint (k = 20 {kcal/mol)/A?). We performed two additional simulations: 

A Starting from the final snapshot of the forward procedure explained above we perform a backward procedure 
with the softer restraint. 

B We repeated both the forward and the backward procedures using the softer restraint 

Structural analysis is shown in Figure 5, where it can be appreciated that only the simulations with protocol A (in 
light blue) were able to correctly form the native WC pairs. Although the restraint stiffness could affect the result, 
we believe that here the differences are mostly due to the stochastic nature of MD. 

Using the snapshots from protocol A, we performed another US simulation. The resulting PMF profiles are shown 
in Figure 6. Also in this case the free-energy landscape is incompatible with the one obtained from the forward 
protocol (compare with Figure 3), indicating convergence issues in the US calculation. 
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FIG. 6: Potential of mean force with alternative backward protocol Potential of mean force (PMF) as a function of 
the distance between the centers of mass of the L2 and L3 loops. Results for Apo form are shown in red (Bwd (k = 20)) 
as obtained from the control umbrella sampling simulations discussed in the Appendix, and aligned with the other two Apo 
profiles described in the Methods section (Fwd in orange, Bwd (k = 40) in black) for comparison. 
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